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1. Introduction 

Numerical solutions of the parabolized Navier- 
Stokes (PNS) equations have been used to obtain a 
better understanding of the qualitative and quanti- 
tative physical phenomena in steady supersonic and 
hypersonic viscous flows. The noniterative numer- 
ical schemes used to solve the PNS equations have 
either required adjustment of solution-dependent 
coefficients for capturing shocks or been inefficient 
on vector supercomputers. The purpose of this 
study is to develop and apply a numerical scheme 
which (1) eliminates the need to make adjustments 
for shock capturing and (2) efficiently utilizes vector 
supercomputers for accurately solving the PNS equa- 
tions for complicated hypersonic flow fields over re- 
alistic vehicle configurations. 

The subject of this study is a finite-difference, 
two-stage, explicit, upwind algorithm for the direct 
(noniterative) integration of the three-dimensional 
PNS equations in a generalized coordinate system. 
The advantages of this type of algorithm are that 

1. The use of upwind flux approximations with 
equation sets containing nonlinear hyperbolic 
conservation laws, such as the pressure and 
convection terms in the PNS equations, allows 
shocks to be numerically captured without 
artificial damping terms which the user must 
adjust. 

2. An explicit integration scheme provides an ex- 
tremely efficient numerical method on vector 
or parallel machines for solving systems of 
equations because the dependent variables can 
be explicitly updated using concurrent ma- 
chine operations. 

The new algorithm uses upwind approximations of 
the numerical fluxes for the pressure and convection 
terms obtained by combining flux-difference split- 
tings (FDS) formed from the solution of an approxi- 
mate Riemann problem (RP). The approximate RP 
is solved by modifying the method developed by Roe 
(1981) for steady supersonic flow of an ideal gas. 
Roe’s method was extended for use with the PNS 
equations expressed in generalized coordinates and 
with Vigneron, Rakich, and Tannehill’s (1978) ap- 
proximation of the stream wise pressure gradient. For 
the three-dimensional PNS equations, both fully up- 
wind and upwind-biased approximations of the pres- 
sure and convection flux derivatives are used. The 
upwind-biased flux approximation is formed by us- 
ing an upwind flux approximation in the direction 
normal to a shock wave and a central flux approxima- 
tion in all other directions. The upwind-biased flux 
approximation eliminates a loss of accuracy in the 
numerical solution experienced in three-dimensional 


flow when upwind flux approximations were used 
in directions tangential to a shock wave whose tan- 
gential velocity was negligible. The upwind fluxes 
are used in a two-stage integration scheme that re- 
duces to MacCormack’s (1969) method when the 
FDS terms are identically zero. 

The PNS equations can be integrated (marched) 
in space using either an iterative, a noniterative, or 
a time relaxation scheme. Time relaxation schemes 
retain the time-dependent terms and use time inte- 
gration methods to obtain a steady state solution 
at a streamwise location before advancing in space. 
A noniterative method is usually preferable over ei- 
ther a time relaxation or an iterative scheme since 
the solution at a given streamwise station is obtained 
directly. A noniterative method is important when 
each integration step’s cost is high. This study con- 
siders only the noniterative schemes applied to the 
PNS equations. 

The conservation law form of the PNS equa- 
tions is usually solved using refinements of the 
finite-difference codes of Schiff and Steger (1979) 
or Vigneron, Rakich, and Tannehill (1978). Both 
codes use a noniterative, implicit, approximate- 
factorization, finite-difference algorithm for integrat- 
ing the thin-layer form of the PNS equations. These 
algorithms are based on numerical schemes devel- 
oped by McDonald and Briley (1975) and Beam and 
Warming (1978). These numerical schemes use cen- 
tral differences to approximate the spatial derivatives 
of the fluxes. These implicit schemes also use recur- 
sive operations which are generally more difficult to 
apply efficiently on vector computers than are ex- 
plicit schemes. Gielda and McRae (1986) took ad- 
vantage of the high vectorizing efficiency of a mod- 
ified form of MacCormack’s (1969) method to solve 
the PNS equations on a Cray 1 supercomputer. They 
achieved total solution times that were competitive 
with existing implicit algorithms for certain classes of 
problems. Conventional central difference schemes 
such as Beam and Warming’s and MacCormack’s 
require manual adjustments of artificial damping 
terms to maintain numerical stability and to elim- 
inate nonphysical oscillations in the numerical solu- 
tions around shock waves. Lawrence, Tannehill, and 
Chaussee (1986, 1987) developed an implicit finite- 
volume scheme for solving the PNS equations which 
used upwind differencing of the convection terms in 
areas of supersonic flow and standard central dif- 
ferencing in subsonic regions. Central differencing 
was used in the subsonic region because numerical 
instabilities occurred when upwind differences were 
used. The disadvantage of this approach is that 
it is difficult to vectorize because of the difference 



switching and the implicit integration. The finite- 
difference upwind method developed in this study 
can be used throughout the flow field and does not 
require any special switching of differences in the sub- 
sonic regime. 

This report presents a new three-dimensional, 
noniterative PNS solver which combines the compu- 
tational speed and second-order marching accuracy 
of a two-stage explicit integration scheme with the 
robust features obtained from upwind approxima- 
tions of the convection terms. The new algorithm for 
solving the PNS equations has the following unique 
features: 

1. Use of upwind approximations of the convec- 
tion terms in the subsonic region 


2. Application of a two-stage integration scheme 
with upwind flux-limited approximations of 
the convection fluxes 

3. A cubic equation defining Vigneron’s splitting 
coefficient in terms of the dependent variables 

4. Use of different upwind flux approximations 
in each stage of the integration algorithm 

The outline of this study is as follows: in section 2 
the PNS equations are derived; in section 3 the inte- 
gration of hyperbolic conservation laws in multistage 
explicit schemes is investigated and applied to two- 
and three-dimensional PNS equations; and in sec- 
tions 4 and 5 solutions using the new algorithm are 
obtained for two- and three-dimensional flows. 
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2. Governing Equations 

The design process for aerospace vehicles has been improved by using computational fluid 
dynamic (CFD) computer codes to better understand the qualitative and quantitative physical 
phenomena in the flow field. Future aerospace vehicles, such as the National Aero-Space Plane 
(NASP), will cruise at hypersonic speeds. At these speeds, the shock waves and boundary 
layers are merged over portions of the vehicle (fig. 2.1) and cannot be solved independently. 
Ideally, one should solve the full Navier-Stokes (NS) equations. However, the memory and 
speed requirements to obtain full NS solutions for a complete vehicle are beyond the capabilities 
of modern supercomputers. The parabolized Navier-Stokes (PNS) equations approximate the 
full NS equations and are more amenable to efficient numerical solution for steady hypersonic 
flow problems. This efficiency arises from spatial rather than temporal integration. Spatial 
integration reduces the problem by one dimension and thus provides a significant savings in 
computer memory and time for a given case. PNS numerical solutions agree with NS numerical 
solutions for steady supersonic and hypersonic viscous flow problems that do not have a strong 
upstream influence from points downstream. In this section, the strong conservation law form 
of the PNS equations is presented for use with a generalized coordinate system. 



2.1. Navier-Stokes Equations 

The flow of a Newtonian fluid can be described by the NS equations, the continuity of mass 
equation, and the energy equation. These equations expressed in Cartesian conservation law 
form are 

Ut + E x + F y 4- G z = 0 (2*1) 

where subscripts x, y, and z indicate partial differentiation with respect to Cartesian coordinates, 
and subscript t with respect to time. The vectors are defined as 

U = [p, pu, pv, pw, e t ] T 
E = Ei - E v F = Fj - F„ G = Gf - G„ 

Ej = [pu, puu + p, puv, puw, ( et + p)u] T 

Et, = [0, T xx , Txy , T xz , UT XX + VT xy + WT XZ - q x ] T > ^ 

F i = [pv, pvu, pvv+p, pvw, ( et+p)v] T 
F y = [0, T X y, 7yy , TyZ) UTxy ~K UTyy “h WTyz Qy\ 

G i — [j pw , pwu , pwv, pww + p , (et +p)w] T 
Gv = [0, T XZ , T yz , T ZZ , UT XZ + VT yz + WT ZZ - q z ] T j 

The fluxes are separated into inviscid (subscript i) and viscous (subscript v) components and p 
is the density; u : v, and w are Cartesian velocity components; et is total energy; p is pressure; 
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r is viscous stress; and q is heat flux. The total energy is defined by 


et= P 



V? + V 2 + w 2 ^ 


(2.3) 


where e is the internal energy. The variables have been nondimensionalized using the following 
relationships: 


t = 


u = 


P = 


L/U 0 

u* 


Un 


PooUl 
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V = 


u oo 

m T * 
T= — 

Too 
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y = T 


w = 
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w 

u r 

e 

w 

u oo 


oo 
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Z = 


P = 




Poo 
Poo ) 


(2.4) 


where T is temperature, \x is viscosity, U is total velocity, L is a characteristic length, * denotes 
dimensional quantity, and subscript oo represents dimensional reference conditions. 


2.2. Generalized Transformation 

The simulation of a steady supersonic or hypersonic flow field with a dominant flow direction 
is the problem of interest here. The streaimwise direction is defined as being aligned with the 
£-axis. The crossflow plane is defined by the coordinates 77 and C- Equations (2T) and (2.2) are 
now changed to a generalized coordinate system using the following transformation: 

£ = £(z) v = v(x,y,z) ( = C(x,y,z) (2.5) 

The indices i, j, k identify discrete points in the £, rj, C computational coordinate system. The 
transformation and the formulas for the metrics are given in appendix A. 

2.3. Parabolized Navier-Stokes Equations 

The PNS equations have evolved from the pioneering work of Rudman and Rubin (1968), who 
derived and numerically solved a composite set of equations valid in both the inviscid and the 
viscous region for steady supersonic and hypersonic flow. Researchers have used various forms 
of composite equation sets to solve steady supersonic viscous flow problems. Common to each 
of the equation sets is that all viscous stresses and heat fluxes in the streamwise direction are 
dropped and the subsonic streamwise pressure gradient is modified or dropped. Each numerical 
scheme obtains a solution at a given streamwise plane, by either time relaxation, iteration, or 
direct solution, before proceeding forward in the streamwise direction to the next plane. 

The PNS equations can be obtained from the transformed NS equations by assuming steady 
flow and by neglecting the streamwise diffusion terms. The transformed PNS equations in strong 
conservation law form are 
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where subscript £, 77, £ indicate partial differentiation with respect to the generalized coordinates 
and J is the Jacobian of the transformation. The streamwise pressure gradient has been split 
between the two vectors E* and P using the technique developed by Vigneron, Rakich, and 
Tannehill (1978). The streamwise inviscid flux is split as 


Ei = E* + P 


(2.7) 


where 

E* = [pu, puu + up, puv , puw, ( et + p)u] T 

P — (0, (1-0 ,)p, 0, 0, 0] T 

The value of Vigneron’s coefficient u varies from 0 to 1. The determination of u is covered in a 
later section. 

The shear stress and heat flux terms after the transformation and the parabolizing assump- 
tions become 


2 n ' 

T ix = [2 (VxUj] + (i«() - ( VyV v + CyV(;) ~ (VzW v + £ 2 t^)] 

T xy = r 7 ~ { r ly u i] + Cy u C + VxVt) + Qx v q) 
aIGqo 


qx (7 - lJM^ReooPr + 

T yy = gj^ [2 (Vy v rj + C y v <;) ~ (VxUy + Cx"U() 

Txz = ( VzUr, + CzU( + VxW v + CxW() 

rieoo 


(vzWj) + G^)] 


> 


Qy - (7 - l)M2,ReooPr ^ + ^ 

T zz = [2(ri z w v + C, z Wq) - (T) x Ur, + £ I«() - (VyV v + £y^)] 

T yz = (TjyWrj + CyW^ + T] Z V V + ( Z V C ) 

rtCoo 

qz = (7 - ljM^ReooPr ^ + ^ 


( 2 . 8 ) 


The Reynolds, Prandtl, and Mach numbers are denoted by Re, Pr, and M. The ratio of specific 
heats is denoted by 7 and the viscosity is calculated using Sutherland’s equation: 


where 



(2-9) 


The perfect gas relationships are used to completely define the system of equations: 

T _ 

P 


(2.10) 


Note that the equations presented here differ from those used by others (Chaussee et al. 1981) 
in that the stress terms are retained in the flux vector E when the outer derivative is other 
than 
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2.4. Geometric Conservation Law 


The second and third lines of equation (2.6) contain the geometric conservation law (GCL) 
(Anderson, Tannehill, and Pletcher 1984) and metrics associated with returning the governing 
equations to strong conservation law form. The sum of the metrics in the GCL equals zero when 
they are analytically evaluated for the full NS equations. In a numerical scheme, the GCL terms 
may or may not be zero depending on how the full NS equations are approximated and how 
the equations and the metrics are differenced. Gielda and McRae (1986) have shown that for 
MacCormack’s (1969) method, the GCL terms are not zero for the PNS equation set (eq. (2.6)) 
for any combination of possible differencing of the metrics. This fact, which is true for all PNS 
solvers, requires that the GCL terms be evaluated numerically as part of the integration scheme 
to cancel nonphysical source terms. 


2.5. Treatment of the Streamwise Pressure Gradient 


In equation (2.6), the streamwise pressure gradient is split between the left and the right 
side of the PNS equation using Vigneron’s coefficient w. If one is interested in using a space 
marching method for integrating the PNS equation set, then the inviscid eigenvalues have to be 
real and the viscous eigenvalues have to be nonnegative and real. These conditions are true for 
supersonic flow. A linear stability analysis by Vigneron, Rakich, and Tannehill (1978) (and later 
extended by Davis, Barnett, and Rakich 1986) shows that the inviscid eigenvalues are real and 
the viscous eigenvalues are real and nonnegative for subsonic flow if a fraction of the streamwise 
pressure gradient is retained. This fraction is obtained by defining u> as 


1 (M { > 1) 



7 Ml 

l+(7-l)M| 


(M* < 1) 


( 2 . 11 ) 


where the axial Mach number is denoted by M £. The coefficient to has to be applied with a 
safety factor a: 

u = min(l, aw) (2.12) 


so that the eigenvalues of the inviscid PNS equation set are real and the equations are hyperbolic. 
The effect of a on w is shown in figure 2.2. 

When u is less than 1, the source term P physically represents the mechanism which would 
allow for upstream propagation of information (Davis, Barnett, and Rakich 1986). Lubard and 
Helliwell (1974) have shown that when P is included in a finite-difference method as part of a 
backward difference, the method is unstable for small marching step sizes. If the source term 
P is dropped, the finite-difference method is stable up to the allowable marching step size for 
the numerical scheme. Therefore, the vector P on the right side of equation (2.6) is dropped 
in numerical calculations. Davis, Barnett, and Rakich (1986) demonstrated that this is a good 
approximation for a high Mach number, weakly interacting flow. However, the source term P 
can be included as a forward difference with a global iteration method (Davis, Barnett, and 
Rakich 1986) for calculating strongly interacting flows after one pass has been made to establish 
an initial pressure distribution. 
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Figure 2.2. Effect of safety factor a on Vigneron’s coefficient u. 

In equation (2.6), cu is included in the differential with respect to £. Once the vector P is 
dropped in the numerical scheme, there is no longer conservation of the streamwise momentum 
(neglecting the streamwise viscous terms) in the subsonic region. Care has to be taken when 
integrating E* so that any variation of u is canceled out in the streamwise (£) direction. If the 
variation of t o is not canceled out, a nonphysical acceleration occurs in the subsonic streamwise 
direction which affects the accuracy of the PNS solution as an approximation to the full NS 
solution . 1 The pressure derivative in the streamwise direction after the vector P has been 
dropped can be expressed as 

= tj-ups + (2.13) 

The last term on the right side has to be canceled out to eliminate the nonphysical accelerations 
caused by the variation of u>. Since u is directly related to the streamwise Mach number, different 
solutions can be obtained for the same problem using slightly different numerical grids if the 
last term is not canceled out, by subtracting its value at the previous solution station from the 
numerical solution. 

2.6. PNS Equations for Use With a Single Pass Method 

The equation solved by a single pass space marching method can be obtained from 
equation (2.6) by dropping the vector P, shifting all the crossflow flux derivatives to the right 

1 The author would like to acknowledge J. H. Morrison, Analytical Services & Materials, Inc. (Hampton, Va.) for 
discovery of this nonphysical acceleration from the variation of u. 
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side, and adding a correction to cancel out any variation of u in the streamwise direction: 


where 


1 

II 

* 

a 

( %E T]yF 

VzG 


GE 

CyF 

C,G\ 


v j 

{ J J 

J 

Jr, \ 

J 

J 

J /c 


+ 


e (t 

V J > 

) + E 

1 

i 1 

3b 

+ F 

c 1-3 
v -h 

+ 

1 1 


+ < 

i 1 

Cj? 

T 

'fell 

"'AJ 

+ JP^ 



(2.14) 


n 

= [o, W, 0, 0 

. or 



(2.15) 


The last term contains a vector which is added to cancel out any variation of u ; in the 
streamwise direction when the source term P is dropped. The PNS equations described by 
equation (2.14) are a mixed set of hyperbolic-parabolic partial differential equations in £-space. 
Given that boundary conditions are known for E* on an 7j-( surface and that appropriate initial 
conditions are known on a surface for £ = 0, the system of equations can be space marched 
(integrated) in the ^-direction. 


2.7. Defining u for Decoding E* 

The primitive flow variables are used in the definition of the fluxes and to display and analyze 
the solution. The primitive flow variables have to be defined in terms of the dependent variable 
E*. Previously, numerical schemes used with the conservative form of the PNS equations have 
required that a change of variables be made from E*, to eliminate the difficulty of decoding 
E* to obtain the primitive flow variables. This difficulty had to do with choosing the sign on 
the square-root function used for determining the streamwise velocity with the steady form of 
Euler’s equations when the flow changes from supersonic to subsonic. Gielda and McRae (1986) 
eliminated this problem by using Vigneron’s coefficient u so that the sign does not change on 
the square-root function when the flow becomes subsonic. They defined the primitive variables 
in terms of E* and u as 


E* = [ pu , puu + up , puv , puw , (e t +p)u\ T = [ E [ , £ £3, £4, E^ T 


(2.16) 


so that 


_El _E\ 

V £* W E{ U 


—b + Vb 2 — 4 ac 
2 a 


,.s 

U 


p = 


£2 - uE{ 

UJ 


> (2.17) 


where 


1 . _ - 7^2 o >(7 - 1 ) 

2 E*[ 2 7 -u,( 7 -1)] [2 7 — w( 7 — 1)] 


B.-l( v z + w A 

L E{ 2 V ^ ) 


The value of u> must be known before the flow variables are computed from equations (2.16) 
and (2.17). Gielda and McRae lagged u in their numerical scheme by defining it as a function 
of the primitive flow variables from the previous decoding of E*. 
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A different approach was taken in this study by defining u> in terms of E* and determining 
it before decoding. This results in a cubic function defining w in terms of E* : 


2(1 + <7)7 2 [C 1 + <7 ) 2 7 2 + j4 2Act7 

+ (7 _ 1)2 w -(7H)5- 


( 2 . 18 ) 


where 


s * — 2 9 

O T7 1 * T7 1 * T?* 

The value of u> can be determined either by solving the cubic exactly or by using Newton- 
Raphson iteration. The advantage of solving for u> from E* is that it allows a larger space 
marching step size in practice, especially when starting from free-stream or approximate initial 
conditions. 
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3. Integration Method 

The objective of this study is to develop a single pass space marching numerical scheme for 
integrating the parabolized Navier-Stokes (PNS) equations which uses upwind approximations 
for the convection terms, is second-order accurate, and can be executed with vector operations. 
In this section the following topics are discussed: an upwind approximation of the convection 
terms using the solution of an approximate Riemann problem (RP), the use of upwind flux 
approximations in multistage explicit integration schemes which have second-order accuracy 
and, finally, a finite-difference upwind algorithm which can be used to integrate equation (2.14). 


3.1. Selecting an Upwind Scheme for the PNS Equations 

Upwind numerical schemes have recently become popular for solving nonlinear, hyperbolic, 
partial differential equations in conservation law form. First-order upwind approximations of the 
flux derivatives are used to eliminate numerical oscillations associated with solutions containing 
discontinuities. This section reviews methods for determining an upwind approximation for the 
purpose of selecting a method for use with the PNS convection terms in an explicit integration 
algorithm. 

An upwind scheme applicable to the steady form of Euler’s equations would be applicable 
to the PNS equations if a modification is made in the subsonic region to handle the splitting of 
the streamwise pressure gradient with Vigneron, Rakich, and Tannehill’s approximation. For 
supersonic flow, the steady form of Euler’s equations is identical to the convective terms of the 
PNS equations. The steady form of Euler’s equations has hyperbolic character for supersonic 
flow and has a more complicated set of eigenvalues and eigenvectors than the unsteady form. 
Most previous applications of upwind schemes have been to the unsteady, conservative form of 
Euler’s equations. 

The characteristics of the PNS convection terms in £-77 space can be examined by first 
determining the high Reynolds number limit (1/Re — ► 0 ) form of equation (2.5) and dropping 
the geometric conservation law (GCL) and source terms: 


E^ + F ir) =0 

E* = — E* 

3 


E*(0, rj) = Eq ( 77 ) 

- _ T] x Ei + T) y Fi + r) z Gi 

- j 


The system of equations is hyperbolic since the Jacobian matrix A is 


(3-1) 



(3.2) 


has real eigenvalues when a < 1 . Equation (3.1) represents an initial- value problem for a 
system of nonlinear hyperbolic conservation laws where the dependent variable E*(£, 7 }) is a 
vector and F is a vector- valued nonlinear function of E . Solutions of equation (3.1) for a 
given set of initial conditions may contain or develop discontinuities of the dependent variable. 
These “weak solutions” of Euler’s equations physically represent shock waves or contact surfaces 
occurring in the flow. Nonphysical weak solutions, called expansion shocks, can also be solutions 
to equation (3.1). The physical solution obeys an entropy condition and is a solution to the 
viscous equations in the limit as e —> 0 : 


E f + 1% - eE 7?7J i e > 0) 


(3.3) 
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The correct physical solution should be predicted by solving the PNS equations since the right 
side of equation (3.3) is approximated by using the physical viscous stresses and heat fluxes. 

The semidiscrete form of equation (3.1) is obtained by approximating the flux on a numerical 
grid in the ^-direction: 


Ej + 


^ F i + i F %- i 


V 


At] 


= 0 


(3.4) 


A first-order upwind approximation of a flux / for a scalar equation at the j + \ point can be 
determined simply by investigating the sign of A, the eigenvalue (or wave speed): 


*i+\ 



(A>0)1 
(A < 0) / 


(3.5) 


For a system of equations where there can be A’s of mixed sign, the upwind flux approximation 
is determined from a splitting of either the flux vectors (FVS) or flux differences (FDS). (See 
summary in Chakravarthy 1987.) Both the FDS and the FVS approach can be used to separate 
the flux into contributions which can be associated with either the positive or the negative A’s. 

The FVS approach is used in the schemes of Steger and Warming (1981) and Van Leer (1982) 
for the unsteady Euler equations. The flux vector is split into two new vectors which have either 
all negative or all positive A’s. A forward or backward difference is applied to the split vectors 
based on the sign of the A’s to obtain an upwind algorithm. Steger and Warming’s method 
is based on a homogeneous property of the Euler equations which does not apply to the PNS 
equations as modified by the splitting of the streamwise pressure gradient. 

The FDS approach is based on information about the evolution of the flow field obtained from 
the localized Riemann solutions between adjacent grid cells. A characteristic decomposition of 
the system of equations defined in the initial-value problem (eq. (3.1)) is made and a Riemann 
problem (RP) is formulated. The initial data for E* are assumed to be piecewise constant 
between points (fig. 3.1). An interface is assumed to exist between the two points at the symbolic 
location of j+ The discontinuous initial data for the initial- value problem define the RP. Since 

E* is a vector and F* is a nonlinear vector, the solution of the RP involves nonlinear algebraic 
equations and logical conditions for determining whether the solution contains a shock wave 
or a smooth expansion. A RP is solved to determine the evolution of the interface in £-space 
and the intermediate values of E*. The solution of the RP (for steady flow of a supersonic gas) 
contains four constant states of E , separated by five waves evolving from the interface (fig. 3.2). 
Each wave is associated with a A of the Jacobian matrix A. The waves can represent a shock 
wave, a rarefaction fan, or a contact surface. Once the intermediate values of E are known, 
the flux difference across a wave can be determined. The flux difference across a wave with a 
positive A is considered a positive flux difference and the flux difference across a wave with a 
negative A is considered a negative flux difference. The positive flux differences are used as part 
of a backward approximation of the flux derivative, and the negative flux differences are used 
as part of a forward approximation of the flux derivative. 


11 



A 



Figure 3.1. Initial data distribution between adjacent points. 

E* 



Figure 3.2. Riemann problem solution. 


Godunov (1960) developed the first upwind scheme for the unsteady, conservative form of 
Euler’s equations. This scheme was based on an exact solution of the RP. The exact solution 
to the RP is expensive to compute since it requires an iteration process. Osher and Solomon 
(1982), Pandolfi (1984), Roe (1981), and others have proposed approximate RP solvers for the 
unsteady Euler equations. Roe (1981) and Pandolfi (1985) have developed schemes for steady 
supersonic flow using the steady form of Euler’s equations. Pandolfi’s schemes are limited 
to Euler’s equations, while Roe’s scheme is applicable to any hyperbolic equation set that 
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has unique eigenvectors. Roe’s method has been applied to the PNS equation set in areas of 
supersonic flow by Lawrence, Tannehill, and Chaussee (1986, 1987). They found that applying 
the upwind scheme in the subsonic region where the streamwise pressure gradient was split was 
“detrimental to both the stability and accuracy of the algorithm.” They used a noniterative 
implicit integration scheme and made a change of variables so that the equations were solved 
in U-space instead of E*-space. For the three-dimensional PNS equations, Lawrence, Chaussee, 
and Tannehill (1987) used the eigenvectors and eigenvalues for Euler’s equations expressed in a 
rotated Cartesian coordinate system for solving the approximate RP in U-space. A modification 
detailed below of Roe’s scheme for use with the three-dimensional PNS equations expressed in 
a generalized coordinate system can be applied for solving the approximate RP in E*-space. 
This modification eliminates the problem experienced by Lawrence, Chaussee, and Tannehill 
(1987) and enables the scheme to be used throughout the flow field, including locations where 
the streamwise. pressure gradient is split. 


3.2. Application of an Approximate Riemann Solver to the PNS Equations 


Roe’s (1981) method is based on solving an approximate RP exactly. The solution of the 
RP is used in a numerical scheme to obtain a splitting of the flux differences in the crossflow 
directions. Consider a discrete grid of r)j points in the 77 -direction. An exact solution is sought 
for the following approximate equation (for the convection terms of the PNS equations) between 
the points j and j + 1 : 


E | + AE; = 0 


(3.6) 


with initial conditions, 


£*( 0 , 77 ) 


%+i (n> j + 5) 


where A is a constant matrix (based on local conditions). The discrete matrix A is formed 
using specialized square-root averaging of the primitive flow variables at points j and j + 1 so 
that conservation properties are maintained. The A variables are formed using a locally constant 
value of the streamwise splitting coefficient u>. The matrix A evaluated at the interface location 
(j -h 5 ) is formed from 



d (^Ei + jFi + G i) 



(3.7) 


The matrix A has the following conservative property if it is evaluated using specialized square- 
root-averaged variables (Roe, 1981): 


F l+i 




(3.8) 


where 






■ 1 + 
J+2 




and where A; is a dummy index for j or j + 1 and the metrics are held constant. For the PNS 
equations, recall that E* has been modified to include only a fraction of the streamwise pressure 
gradient in the subsonic region. The allowable amount of the streamwise pressure gradient can 
change rapidly in the 77 -direction. This variation of d was eliminated in the projection of E* 
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into F in the definition of the RP between points j and j + 1 by defining E* using a locally 
constant value of u: 



n 

pu, pu +UJ- + \P, puv, puw, ( et + p)u 


(3.9) 


This modification was necessary to eliminate the problem experienced by Lawrence, Tann ehill, 
and Chaussee (1986, 1987) in the subsonic region. The requirement of using a fixed value of u> 
for the RP increases the difficulty of using an implicit integration scheme when u> ± 1 (subsonic 
region). This is because the E* at point j in the RP defined at j + \ is not equal to the E* at 
point j defined in the RP at point j — \ because of the different values of ui used at the j ± \ 
points. 

The objective of solving the RP between points j and j + 1 is to split the flux difference in 
equation (3.8) into five parts, one for each eigenvalue, which can be associated with either the 
positive or the negative eigenvalues. Roe’s method for solving the RP consists of first calculating 
the square-root-averaged variables for the interface at j + The eigenvalues A and eigenvectors 
e of the Jacobian matrix A, and the wave strengths a are calculated with the square-root- 
averaged variables using the equations given in appendix B. The wave strengths are defined 
as 

£ KA»> J+ i = (7) . , (%Vi - %*) (3.io) 

m= 1 ' w + S 

The flux difference across the rath wave is A m a m e m . The sum of the flux difference across all 
the waves is equal to the difference of the flux between points j and j - hi. 


5 

T. l = Fj+1 ~ Fj 

m— 1 


(3.11) 


The total flux difference between points j and j + 1 can now be split into the total positive (df + ) 
and negative (df - ) flux differences. The df + and df - vectors are calculated from 


? i+ i - Fj = dft , + df - , 
] J+h J+\ 


(3.12) 


where 


df + j = y 

rn—1 


\n + |A m| 


( a m&m)j + l 



EA m |A m | , A . 

2 

m= 1 ^ 


(3.13) 


The eigenvalues and eigenvectors for the three-dimensional inviscid PNS equations in generalized 
coordinates were determined in part using the symbolic manipulation language MACSYMA. The 
vectors df and df are used as the building blocks for obtaining an upwind flux approximation 
at the j + j point. 


3.3. Upwind Flux Approximations 


A first-order upwind flux approximation at j + \ can be obtained by modifying either a 
forward, a backward, or a central flux approximation with FDS determined from the solution 
of the RP: 







(3.14) 
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A numerical scheme that uses a first-order upwind approximation of the flux has the advantage 
of resolving discontinuities without spurious oscillations. Unfortunately, the dissipation inherent 
in the first-order upwind scheme makes it impractical for global use. In practice, second- and 
higher-order flux approximations are used with numerical schemes to minimize truncation errors. 
The higher-order flux approximation cannot resolve a discontinuity such as a shock wave without 
an overshoot or undershoot. The oscillations around a shock wave are minimized either by 
adding additional dissipation to the numerical scheme or by modifying the flux approximation 
with a flux limiter. 

The classical way of minimizing the oscillations is to use additional dissipation. Dissipation 
is added to the numerical scheme by including either a second- or a fourth-order derivative of the 
dependent variable multiplied by a user-specified constant. The disadvantages of this procedure 
are determining the best value for the constant and adding a nonphysical stresslike term to the 
equation set that is being solved. 

Another approach is to use a nonlinear method to change or “limit” the higher-order flux 
approximation to first-order in the neighborhood of a discontinuity to eliminate numerical 
oscillations (fig. 3.3). A second-order upwind approximation of the flux using flux limiters 
is 



+ 5 


(<> - df- + ,) 


(3.15) 


The flux differences designated with an overbar are treated with the minmod flux limiter. 

df+ i = minmod ( df + i , 

3 2 V j-j 

df . , 3 = minmod ( df“ 

J+v \ J+i 

The minmod flux limiter is defined as 



minmod (x,y) = sign(x) max {0, min [|x|, y sign(x)]} 


(3.17) 


The minmod function limits the overbar flux differences to a first-order approximation around 
captured shocks in order to minimize oscillations. The parameter (3 varies slightly for different 
flux approximations. For second-order upwind flux approximations, (3 has a maximum value 
of 2 for obtaining oscillation-free shock capturing. The advantage of the nonlinear flux limiter 
is that oscillation-free results can be obtained without adding artificial stresses to the numerical 
scheme. Different types of flux limiters are in use and the reader is referred to the work of Sweby 
(1984) for more information. 
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Figure 3.3. Comparison of unlimited and flux-limited shock capturing with a second-order flux approximation. 

3.4. Second-Order Explicit Upwind Integration Method 

First- and second-order upwind flux approximations for the PNS inviscid equations in 
£-7/ space were defined in the previous sections. This section presents the integration of the 
initial- value problem described by equation (3.1) using a method that incorporates upwind flux 
approximations and is globally second-order accurate. The flux approximations in MacCor- 
mack’s (1969) scheme are modified to obtain a flux-limited version of Warming and Beam’s 
(1975) second-order explicit upwind scheme. 

MacCormack’s (1969) method is a two-stage explicit scheme that has been used extensively 
for solving the Euler, full Navier- Stokes, ‘and, recently by Gielda and McRae (1986), the PNS 
equations. The one-sided inviscid flux approximations used in MacCormack’s method are 
modified with the FDS obtained from the solution of the RP. The resulting unlimited form 
of the algorithm is similar to the Warming and Beam upwind (WBU) algorithm. The WBU 
scheme has twice the linear stability limit of MacCormack’s. The larger integration step size 
of the WBU scheme compensates for some of the additional cost of determining the solution to 
the RP’s. The MacCormack and WBU schemes are classified as Lax and Wendroff (1960) type 
schemes. 

Second-order and higher integration schemes for solving initial-value, boundary-value prob- 
lems are derived using either the fully discrete method of Lax and Wendroff (1960) or the 
semidiscrete method of lines. The Lax- Wendroff scheme is a finite- difference method that is 
derived by satisfying a Taylor series expansion about the solution point. The derivatives with 
respect to the independent variable of expansion are replaced with the original partial differen- 
tial equation. Another way to integrate the equation set is the method-of-lines approximation. 
The first step in the method-of-lines procedure is to express the partial differential equation 
in semidiscrete form. The remaining partial derivatives are treated numerically as ordinary 
differential equations. Two- and three-dimensional problems with source terms can be easily 
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handled, and any appropriate ordinary differential equation solver can be used to perform the 
integration. 

The method-of-lines procedure can be modified for determining a numerical scheme equiva- 
lent to the MacCormack and WBU schemes. First apply a second-order modified Euler method, 
also called Heun’s method (Gear 1971), to the initial- value problem expressed in semidiscrete 
form (eq. (3.4)): 




= E f - 


Jf - 

F {0) 


a?7 \ 
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(3.18a) 

(3.18b) 


where the superscript represents the stage of evaluation and (0) represents the initial value. The 
method-of-lines procedure guarantees a global second-order scheme when the midpoint fluxes 
are approximated with the same second-order flux approximations in each stage. Equivalent 
Lax-Wendroff methods can be formulated if different types of flux approximations are used in the 
different stages. The accuracy of the method must be checked when different flux approximations 
are used unless the scheme is a known Lax-Wendroff method. 

MacCormack’s scheme is obtained using one-sided flux approximations: 


Stagel F i°Ji =F i 0) 

Stage 2 Ff =F; 0) 

j +2 J 


f { i; ) = F.^ 

3+y 


(3.19) 


l j+ 1 


Substituting for the flux differencing evaluated at level (0) in the second stage (eq. (3.18b)) 
with the first stage equation (eq. (3.18a)) results in the traditional form of the second stage of 
MacCormack’s algorithm: 


E ll) = B® - (f! 0) - f! 0) ) 

3 3 A 77 V *i- 1 / 


1 

fe <0, + Ef-£ 

fFj 1} - fJ 1} ) 


2 

3 3 At/ 

V l j + 1 l 3 ) 

\ ) 


The upwind scheme is obtained by adding FDS to the one-sided flux approximations: 


Stage 1 
Stage 2 


p(°) 

pK°) 

l j+\ 


= Fj 0) + df(° 1 ) 

= f„ ( 0) + df :f? + 


3+1 


+( 0 ) 
1 

1 


■Ki) _«Ki) 


l i+l 


V ) 1 -df 


+( 1 ) 


l j+ 1 


3+1 


1 (m; l ‘ 

2 \ 3- 

1 f^ +{0) 


^ 1 ) 


(3.20) 


(3.21) 


The flux approximations used in stage 2 differ from those in stage 1 and contain terms evaluated 
at both levels. The numerical fluxes approximate the flux at a particular point in space, which 
may differ from stage to stage. The flux definitions in stage 2 are not unique. For example, 
all the terms evaluated at (0) could be placed together. The above form was preferred since 
it represented a second-order discrete approximation of the flux in 77 -space. Different flux 
approximations used in different stages of Runge-Kutta integration schemes have been used 
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in developing the third-order schemes of Rusanov (1970) and Burstein and Mirin (1970). The 
traditional form of the second stage of the WBU scheme is obtained by substituting the equation 
for the first stage into the second: 


e* ( 1) = r (0) - — 

J J At? [ 


F- 0) - fJ 0) + ( df 

l J l J~l 


5 + e ; (1) - 1 { f s. - - <y ) 


(3.22) 


The advantage of this upwind scheme is that it can easily be obtained by modifying the inviscid 
fluxes used in MacCormack’s scheme. The above scheme differs in practice from other upwind 
schemes because the first-order upwind flux approximation is based on modifying either a 
forward or a backward flux instead of a central flux (eq. (3.14)). The metrics used in the 
definition of the FDS are defined so that any downwind contributions exactly cancel. 

Consider the first-order approximation of the flux derivative including the geometric conser- 
vation law (GCL) terms on a grid where if = rj(y) with the upwind flux definition based on a 
backward flux approximation: 



If all the eigenvalues are positive, the FDS terms would be zero and result in a backward 
difference approximation of the flux derivative. If all the eigenvalues are negative, the FDS 
terms exactly cancel out the backward difference if 



(A < 0) 


(3.24) 


and thus 

dt J-5 = (?),_■ ( F <* - F <»-.) < AU A < °> 

which results in a forward difference approximation of tide flux derivative. 

If an upwind flux is formed based on a forward flux approximation, then the forward 
difference would be canceled by the FDS terms if all the eigenvalues are positive and the metrics 
are defined by 

(y), + . • (tW ; (A > 0) 11,251 

and thus b 

d £i = 0O i+ i - (AUA>0) 

The metrics for a general transformation, i ? = t](x, y, z), would be defined at the same points 
as above for the positive and negative flux differences. The different approximations for the 
upwind flux should be alternated after one complete cycle of the algorithm to eliminate biasing. 
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3.5. Viscous Stress and Heat Flux Differencing 

The finite-difference approximations of the stress tensor and the heat fluxes are made so that 
their derivatives formed in each stage of the algorithm are second-order accurate. This requires 
separate approximations of the stress tensor and the heat fluxes for each differencing direction. 
The separate approximations of the stress tensor and heat fluxes result in a source term that is 
canceled by including the GCL terms. The derivatives of the stress tensor for the x-momentum 
equation including the first GCL term are 

^hjxx Vy T xy flzTxz ^ Cy T zy Cz T xz ^ .// ^Cx 

where two discrete approximations of the same stress are used. The single prime denotes stresses 
to be differenced in the 77 -direction, and the double prime the C-direction. Let j represent the grid 
points in the 77 -direction, and k the C-direction. For example, assume that a forward difference 
of the viscous stresses in the first stage of the algorithm is required. The discrete approximation 
of T xy in a generalized coordinate system (A 77 = AC = 1) would be 


- . . . (3.26) 

C 



(3.27) 


This differencing scheme for the stress and heat flux terms results in a source term that has 
to be canceled by including the geometric conservation law (GCL) in the integration algorithm 
(see Gielda and McRae 1986). 

3.6. Explicit Upwind Integration Scheme for the Three-Dimensional PNS 
Equations 


An explicit upwind integration scheme for the three-dimensional PNS equations can be 
formed using the flux-limited form of the WBU scheme and the stress and heat flux differencing 
developed by Gielda and McRae (1986). A second-order, two-stage, explicit, upwind scheme for 
the PNS equations is 
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The superscript n stands for the values determined at the initial condition, the p values are 
based on the result of the first stage, and n + 1 stands for the solution point. The fluxes F„ and 
G v denote the diffusion terms. The single prime denotes viscous stress and heat fluxes to be 
differenced in the 77 -direction, and the double prime denotes differencing in the (-direction. 


3.7. Geometric Conservation Law 


The GCL term included in the numerical algorithm is used to cancel source terms occurring 
from the differencing of the stress terms as mentioned previously and the use of the Vigneron 
coefficient in the dependent flux vector (see Gielda and McRae 1986 for details). The GCL term 
in equation (3.28) is defined as 



(3.29) 


3.8. Implementing Boundary Conditions 

Boundary conditions at the wall are implemented by applying the no-slip condition at the wall 
to the velocity components, applying the constant wall temperature condition, and satisfying a 
normal momentum equation solved for pressure at the wall. The density at the wall is determined 
from the wall temperature, pressure, and equation of state. The discrete equation for the wall 
pressure was established by applying the no-slip velocity condition to the y- and z-momentum 
equations where 77 is the direction that intersects the wall. The y- and z-momentum equations 


20 



evaluated at the wall are 


VyPr) + CyPc — Vx(^xy)r} + Vy( r yy)rf + Pz{^zy)r\ 

+ C x(jxy)(^ + (y( r yy)(; + Cz( r zy)^ 
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The derivative is eliminated between the two momentum equations: 
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(3.31) 


where 


V = Tf x U + TjyV + 7] Z W , IT = Cx^ + Ct/^ + Cz™, 1? = Cz^ - Cy™ 


A second-order, finite-difference module for equation (3.31) is solved for the wall pressure. The 
above formulation differs from the more approximate method obtained by neglecting the stress 
derivatives. Including the stress terms in equation (3.31) eliminated severe oscillations of the 
surface pressure experienced at the tip of delta wings in high Mach number flows. 


3.9. Selecting a Marching Step Size 

The allowable step size is calculated from a combination of the maximum allowable inviscid 
and viscous step sizes for a linear system. The inviscid upwind algorithm has a Courant- 
Friedrichs-Lewy (CFL) linear-stability limit of 2, and a viscous limit of V2 for two-dimensional 
problems and Va for three-dimensional problems. If the distance between points in the 17- and 
C-directions is denoted as A r and As, respectively, a marching step size Ax can be calculated 
by neglecting the crossflow velocities and assuming that the first point off the wall is subsonic: 

/\ d / \ 

S M 1 , 2 ". ( ’ 

2min (Ar, -c> + l"‘ Alt 


where 


A d = 


Ar As 
\/Ar 2 4- As 2 


Equation (3.32) is usually applied with a safety factor of about 0.90 for simple geometric shapes. 
Note that this equation gives a marching step size that is twice that of MacCormack’s scheme. 
The starting step size should be calculated by replacing the local axial Mach number M ^ with 
one slightly larger than the desired M ^ for one point off the surface. 


3.10. Three-Dimensional, Upwind-Biased Method 

The PNS equations are integrated using space marching procedures. This allows a two- 
dimensional flow problem to be solved with a one-dimensional numerical method, and a 
three-dimensional flow problem with a two-dimensional numerical method. The upwind 
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flux approximations for the pressure and convection terms are formed from one-dimensional 
numerical operators. The one-dimensional numerical operators are determined analytically by 
solving a RP. The Riemann solution is for an approximate form of the two-dimensional inviscid 
PNS equations which locally models the inviscid flow between two points. Unfortunately, there 
is not a two-dimensional analytical solution of steady three-dimensional flow which can be 
used in determining an appropriate two-dimensional numerical operator for approximating the 
numerical fluxes. In practice, the numerical fluxes for a multidimensional numerical problem 
are usually determined by a sequence of one-dimensional numerical operators. Algorithms for 
the three-dimensional PNS equations require two different flux approximations in the crossflow 
plane. The fully upwind algorithm for the three-dimensional PNS equations presented earlier is 
based on two one-dimensional upwind operators determined from the Riemann solution for two- 
dimensional steady flow. The combination of two sets of RPs, one for each coordinate direction 
in the crossflow plane, does not always correctly model the physics of a three-dimensional flow 
field. The purpose of this section is to investigate alternate ways of approximating the convection 
fluxes in the crossflow plane to model the physics of the flow field more accurately. 

Upwind flux approximations constructed from analytical solutions of RP’s are used because 
of their excellent shock-capturing capability. For each integration step of a one-dimensional 
upwind numerical method, a shock wave is captured between only two points. This represents 
a point discontinuity along a line of data. For three-dimensional steady flow, the shock wave 
represents a discontinuous curve in the crossflow plane. The use of one-dimensional analytical 
solutions to determine upwind flux approximations for multidimensional schemes can result in 
difficulties associated with the discretization of a continuous problem containing a shock wave. 
When the shock wave is mapped onto a discrete grid for the crossflow plane (fig. 3.4), the 
discontinuity cuts across data cells unless the shock aligns with the numerical grid. The fully 
upwind algorithm treats the discontinuous data as jumps and attempts to model a shock at each 
cell interface in all directions. In this case, the upwind flux approximation does not represent 
the physics as clearly as in the two-dimensional case. When the crossflow velocity tangent to 
a shock wave is near zero, the numerical dissipation of the scheme in this direction is small. 
The numerical dissipation that is inherent in the numerical method depends on the magnitude 
of the eigenvalues. Three out of the five eigenvalues are directly related to the magnitude 
of the crossflow velocities. Significant errors can occur when the velocity tangent to a shock 
is near zero and a slight difference occurs between adjacent cells next to the captured shock 
(fig. 3.5). The numerical results may drift away from the correct solution when a shock wave 
is not perfectly aligned with the numerical grid in areas where the numerical scheme has little 
natural numerical dissipation. As a shock wave moves across a numerical grid, errors are made 
as the discontinuity jumps across finite areas. If the velocity tangent to a strong shock is nearly 
zero, small disturbances can occur in the tangential velocity that transports large amounts 
of conserved quantities (Colella and Woodward 1984). These errors are minimized in regions 
where there is a large enough tangential velocity providing enough dissipation from the numerical 
method to ensure the correct solution. This problem has been handled by adding dissipation 
to the numerical method using a number of different devices. A commonly used procedure 
is a device originally proposed by Harten (1983) for eliminating expansion shocks. Harten’s 
device maintains the absolute value of the eigenvalues above a minimum value to guarantee a 
minimum amount of numerical dissipation in all flow regions. Specifying the correct amount is 
difficult when the viscous terms are present since the convection eigenvalues are near zero in 
the boundary layer. Colella and Woodward (1984) developed a number of different methods 
for eliminating the multidimensional oscillations with their RP solver for the Euler equations. 
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Unfortunately, each one of these methods uses coefficients that require tuning by the user. The 
main reason for using an upwind scheme was to eliminate the need to use additional dissipation 
that would depend on user-specified parameters. 




Shock interface 


Uneven shock 
capturing 


Figure 3.5. Uneven shock capturing in the crossflow plane. 


An alternate approach to adding dissipation explicitly is to use a numerical method that 
models the physics of the flow more accurately. A two-dimensional numerical scheme is defined 
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as a scheme which selects different flux approximations based on some physical characteristic 
of the flow field that does not occur in the one-dimensional scheme. Murman and Cole (1971) 
and Jameson (1974) used two-dimensional schemes for integrating the potential equation in 
regions of transonic flow. Murman introduced the idea of equation-type-dependent differencing 
with the use of upwind differencing in hyperbolic regions and central differencing in elliptic 
regions. Jameson (1974) applied a rotated differencing scheme to the potential equation that 
used upwind differencing only in the local stream direction when the equation set is hyperbolic 
(for supersonic flow) and central differencing elsewhere. Davis (1984) developed a rotational 
upwind-biased difference scheme for use with Euler’s equations. The first step of Davis’ scheme 
was to determine the direction in which a shock wave (i.e., RP) would be oriented if it existed 
in the finite-difference stencil (fig. 3.6). An upwind algorithm would then be used in a local 
coordinate system rotated to be normal to the orientation of the assumed shock wave. Roe (1986) 
proposed a two-dimensional RP algorithm which has not been implemented. Hirsch, Lacor, and 
Deconinck (1987) have investigated a method in which the two-dimensional unsteady form of 
Euler’s equations can be diagonalized for special flow conditions. 



Figure 3.6. Shock angle in the computational plane. 


C 


The two-dimensional problem described previously has to do with the mapping of a shock 
wave as a discontinuous curve. Riemann methods are based on the assumption that the data are 
piecewise continuous. This assumption is correct in the direction normal to the shock wave. The 
data in the tangential direction to the shock wave should be continuous without jumps. Davis’ 
(1984) scheme is based on this observation. Recall that the fully upwind scheme described in 
section 3.6 becomes a central difference scheme based on MacCormack’s (1969) method when 
there is no jump in the initial data used in solving the RP’s. A combination method based on 
Davis’ ideas can be incorporated by using the solution of the RP only in the directions normal 
to an assumed shock wave position. In this case the flux approximation at an interface now 
depends on the angle of an assumed shock wave. Since the RP is solved for an interface in the 
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computational plane, the angle that the assumed shock wave would make at this interface has to 
be determined. The assumed shock angle can be determined using geometry by recognizing that 
there is no velocity jump in the direction tangential to a shock wave (Davis 1984). A quadratic 
equation for the shock angle derived from the geometry in the computational plane is 

F c tan 2 $ - tan 4* (TF C -V v )-W v = 0 (3.33) 


where V and W are defined as for equation (3.31). Equation (3.33) differs from Davis’ 
formulation of the shock angle since it is not based on an approximation. Once the angle 
of the shock in the computational plane at the interface is known, then the FDS terms can be 
modified by 
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The parameter 6 can vary from 0 to 1. For the interface location between points in the 77 - 
direction, 6 = |cos$|. For the interface location between points in the C-direction, 6 = |sin$|. 

To summarize, first the RP’s are solved in the computational coordinate system as if fully 
upwind differencing is to be used for the pressure and convection terms. Second the angle at 
which a shock wave would occur at the interface is calculated from equation (3.33). Finally the 
FDS expression is multiplied by either the sine or the cosine of the shock angle for the interface 
(eq. (3.34)). Following this procedure, an upwind differencing stencil is used for the pressure 
and convection terms in the direction normal to a shock wave and MacCormack’s differencing 
stencil is used in the direction tangential to a shock wave. 
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4. Two-Dimensional Results 

The explicit upwind algorithm (section 3.6) is 
used to solve the parabolized Navier-Stokes (PNS) 
equations for three two-dimensional supersonic and 
hypersonic viscous flows. The explicit upwind PNS 
solver applied to two-dimensional flow problems has 
the advantage that the numerical problem is one- 
dimensional. The Riemann problem (RP) used in 
this study models the convection terms of the two- 
dimensional PNS equations. The successful simula- 
tion of the two-dimensional test cases validates the 
upwind flux approximations determined from the so- 
lution of the RP. Solutions are calculated for laminar 
flow over a flat plate, across a ramp, and through an 
inlet. The calculations were done on both a Gould 
minicomputer and the Numerical Aerodynamic Sim- 
ulation (NAS) Cray 2 supercomputer. The calcula- 
tions were repeated on the Cray 2 to obtain execution 
times for the vectorized code. 

The numerical grids were generated at each new 
station using algebraic stretching functions to cluster 
the spacing of the points near the wall. The stretch- 
ing transformation can be found in the text of An- 
derson, Tannehill, and Pletcher (1984): 


r(j) = r( 1) + [r(nj) - r(l)] h 



(4.1) 


where 


h(*) = l-0 + 
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The distance away from the wall (j = 1) is denoted 
by r, and rij denotes the total number of points. This 
transformation clusters points around j — 1 when the 
values of the stretching parameter (3 are close to 1. 


4-1. Supersonic Laminar Flow Over a Flat 
Plate 

The first case is for supersonic laminar flow over 
a flat plate. The conditions are 

Moo = 2.0 Re £ = 1.65 x 10 6 L = 1.0 m 

Pr = 0.72 Too = T w = 222 K a = 0.80 

where T w is the wall temperature and a is the safety 
factor introduced in equation (2.12). The supersonic 
laminar flow field for a flat plate at Mach 2.0 is com- 
prised of a weak leading edge shock and a Blasius 
boundary layer profile (fig. 4.1). The correct predic- 
tion of both the velocity and the temperature profile 
verifies that the numerical scheme accurately cap- 
tures the weak shock and resolves the viscous and 
heat fluxes. The use of an explicit scheme for this 


case requires a large number of marching steps be- 
cause of the small normal grid spacing required to 
represent the Mach 2.0 boundary layer. The results 
of the upwind scheme at the final exit station show 
that the numerical dissipation inherent in the upwind 
differencing of the convection term does not adversely 
affect the resolution of the viscous terms. 

The computational domain was sized to include 
the leading edge shock. The vertical height of the 
grid was initially set to y = 0.06 and was made to 
grow linearly to a maximum height of y = 0.525 
at x = 0.93. The computational grid consisted of 
104 points in the 77-plane. The code was started 
using free-stream initial plane conditions with an 
initial maximum marching step size of Ax = 0.00015 
until x = 0.025. Once this point had been reached, 
the step size was varied based upon the result of 
equation (3.32). The grid stretching parameter was 
initially set to 1.01 until after startup was completed 
and then varied to maintain the first point off the 
surface at approximately M = 0.15. The number 
of points for the 77-plane was determined to ensure 
having at least 10 points in the axial boundary layer 
profile at x = 0.93. 

The results of the explicit upwind code are com- 
pared at x = 0.93 with the results of Lawrence, Tan- 
nehill, and Chaussee (1986) who used an implicit 
PNS code (Beam and Warming 1978). The tempera- 
ture and axial velocity profiles at x = 0.93 are shown 
in figures 4.2 and 4.3. The Beam- Warming code was 
started using the results from a boundary layer code 
at x = 0.3 and used approximately twice the number 
of equally spaced points in the boundary layer (only 
half are shown) . The boundary layer profiles compare 
favorably even with the different initial condition and 
stretched grid used with the explicit upwind code. 

This stretched grid was scaled to the boundary 
layer growth to maintain the same Mach number at 
the first point off the wall. In so doing, the scaled 
grid minimized the effect of the source term in 
equation (2.14). To illustrate the effect of two 
additional runs were made on a second grid with 
a fixed point spacing with the same normal point 
distribution as at the final station of the original 
scaled grid. On the second fixed grid, the source term 

is significant because of the streamwise change 
of each grid point’s Mach number as the boundary 
layer grows. Figure 4.4 shows the effect on the 
temperature profile of including and neglecting 
in the numerical scheme when obtaining a solution 
on the two fixed grids. The error of not including 

is clearly seen in the drop in the maximum 
temperature. 

The subroutine used to calculate the flux- 
difference splitting (FDS) required 49 percent of the 
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central processing unit (CPU) time. The algorithm 
required 0.216 x 10“ 4 sec per longitudinal step per 
normal grid point on the Cray 2 computer. 

4.2. Hypersonic Laminar Flow Over a 

Compression Ramp 

The second test case is the simulation of two- 
dimensional hypersonic laminar flow across a flat 
plate and over a 15° ramp. The conditions are 


are defined as 
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Moo = 14.1 Re £ = 1.04 x 10 5 L = 0.439 m 
Too = 72. 2K T w = 297.0 K Pr = 0.72 
a = 0.75 

where the characteristic length, on which the 
Reynolds number is based, is the length of the flat 
plate. The 15° ramp is one of many ramps studied 
experimentally at the above free-stream conditions 
by Holden and Moselle (1970). The distinguishing 
features of this case are the interaction between the 
viscous and inviscid flows and the lack of evidence 
of any separation at the base of the ramp in the ex- 
perimental data. The high Mach number flow on 
the flat plate results in the pressure profile normal 
to the surface being curved in the boundary layer. 
The inviscid flow structure consists of a leading edge 
shock which intersects the induced shock formed by 
the ramp. The two right-running shocks intersect to 
form a single stronger shock, an expansion fan, and 
a contact surface (fig. 4.5). 

The computational grid consisted of 45 points in 
the 77 -direction. The top of the grid was defined by 

_ f 0.135 + 0.104a; (For x < 1.65) 

Vtop | 0.3066 + tan(15°)x (For x > 1.65) 

The grid was stretched in the normal direction using 
a constant stretching factor of 1.04. The law-of-the- 
wall coordinate (y 4- ) for the first point away from the 
wall varied between 0.87 and 0.33 from the leading 
edge of the flat plate to the corner of the ramp and 
reached a maximum of 2.7 on the ramp. The code 
was started using free-stream initial conditions and a 
step size of 0.0005. The step size was reduced when 
the result of equation (3.32) was smaller than 0.0005. 

The results of the explicit upwind code are com- 
pared with the experimental results of Holden and 
Moselle (1970) and the numerical results of Hung 
and MacCormack (1976) for the Navier-Stokes (NS) 
equations. The results for surface pressure (cp), heat 
transfer (c/j), and skin friction ( cj ) coefficients are 
shown in figures 4.6, 4.7, and 4.8. The coefficients 


where 
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and subscript w indicates values at the wall, n repre- 
sents the surface normal direction, and 9 represents 
the angle which the 77 -coordinate makes with the sur- 
face normal direction. A slight overprediction of the 
experimental data can be observed that is consistent 
with numerical results using other methods (Gielda 
and McRae 1986; Lawrence, Tannehill, and Chaussee 
1986). Lawrence et al. attributed this difference to a 
flow misalignment with the ramp. The explicit up- 
wind results agree favorably with the Navier-Stokes 
solution except around the beginning of the ramp. 
Since the PNS equations do not allow upstream prop- 
agation of information in the subsonic region, the 
Cp, c/i, and cj results do not increase until the cor- 
ner of the wedge is reached, while the experimental 
and NS results definitely increase before the corner. 
The comparison of cj in figure 4.8 shows that the 
flow was very close to separating at the corner of the 
ramp. The excellent agreement of the explicit up- 
wind Cf with the Navier-Stokes results before and 
after the ramp corner was attributed to the addition 
of the source term Recall that was added to 
cancel the variation of u in the streamwise direction 
(eq. (2.14)). 

The interaction of the leading edge shock and the 
wedge shock is shown in the Mach number contour 
(fig. 4.9) and the pressure contour (fig. 4.10) from 
the explicit upwind algorithm. The intersections of 
the leading edge shock and the wedge shock, the 
new combination shock, and the expansion fan are 
sharply defined in the pressure contour. The contact 
surface can be seen immediately after the shock in- 
tersection in the Mach contour. Notice that the pres- 
sure contour is oscillation free near the wall. This is 
in contrast to the implicit application by Lawrence, 
Tannehill, and Chaussee (1986), which showed pres- 
sure oscillations near the wall. They attributed the 
oscillations to the switch from their upwind scheme 
in the supersonic region to a central difference al- 
gorithm in the subsonic region. When the above 
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switch was implemented in the explicit upwind al- 
gorithm, the stability and accuracy of the scheme 
deteriorated dramatically. Therefore, it is important 
to use the same differencing scheme throughout the 
subsonic region. 

4.3. Hypersonic Laminar Flow Through 
an Inlet 

The final two-dimensional test case is for a com- 
plicated hypersonic flow field in a converging inlet. 
The geometry of the inlet is given in figure 4.11. The 
free-stream conditions are 

Moo = 15 Re L = 8.0 x 10 4 L = 0.4 m 

= 100 K T w = 1000 K Pr = 0.72 

a = 0.75 


The characteristic length is the length of the flat plate 
at the entrance of the inlet. The NS results (New- 
some, Walters, and Thomas 1987) for this inlet flow 
field indicate a strongly interacting flow which expe- 
riences streamwise separation at the compression and 
expansion corners and at the intersection of the shock 
waves with the wall. One-sweep PNS solvers can 
march through streamwise separation points with- 
out predicting separation. Therefore, a one-sweep 
PNS solution does not accurately represent the NS 
solution in areas around the separation regions. The 
one-sweep PNS solution can be used as initial con- 
ditions for NS solvers or for preliminary evaluation 
of the inviscid flow field structure. The purpose 
of this test case is to evaluate the performance of 
the upwind scheme for an inviscid flow containing a 
number of intersecting and reflecting shock waves. 
This is a numerical test case solved by Lawrence, 
Tannehill, and.Chaussee (1986) and Newsome, Wal- 
ters, and Thomas (1987). Lawrence et al. solved-the 
PNS equations using an implicit noniterative upwind 
scheme. Newsome et al. solved both the PNS equa- 
tions and the thin-layer Navier-Stokes equations us- 
ing an implicit upwind relaxation method. 

The numerical domain for the explicit upwind 
scheme extends in the vertical direction from the 
wall to the centerline of the inlet. The numerical 
grid consisted of 45 points in the vertical direction 
which included the 2 points required for the sym- 
metry boundary condition at the inlet centerline. A 
constant stretching coefficient of 1.08 was used for 
the entire length of the inlet. The normal point dis- 


tribution is identical to what was used by Newsome 
et al. The centerline and the wall boundaries are 


2/top — 0.375 

f 0 


2/wall — \ 


{x — 1) tan(15°) 
[ tan(15°) 


(x < 1) 

(1 < x < 2) 

(x > 2) 


The code was started using free-stream initial condi- 
tions and a streamwise step size of 0.0002. The step 
size was reduced if the result of equation (3.32) was 
less than 0.0002. A safety factor of 0.5 was used in 
the step size calculation to maintain stability when 
the shock waves intersected the inlet walls. This cal- 
culation required 20 220 marching steps to reach the 
end of the inlet (x = 4), using 23 seconds of CPU 
time on the NAS Cray 2 computer. 

Pressure, heat transfer, and skin friction coef- 
ficients from the explicit upwind PNS solution are 
compared with Newsome et al.’s PNS and NS solu- 
tions. The pressure coefficients (fig. 4.12) and the 
heat transfer coefficients (fig. 4.13) agree favorably 
with their PNS solution and NS solution except in 
the regions around the flow separation. A slight dif- 
ference is observed near the location of the second 
shock reflection in the inlet for both the pressure 
and the heat transfer coefficient. This difference is 
attributed in part to the difference in modeling the 
symmetry condition with the finite-difference formu- 
lation used in this study and with the finite-volume 
formulation used by Newsome et al. The compar- 
ison of skin friction coefficients (fig. 4.14) is favor- 
able only for the flat plate region. The complicated 
shock wave structure for the inlet can be observed in 
the pressure contour (fig. 4.15) and the Mach num- 
ber contour (fig. 4.16). The shock waves are sharply 
captured using the upwind scheme without pre- and 
post-shock oscillations. While the overall trends of 
the two PNS solutions are similar, part of the differ- 
ence may be attributed to the more complete mod- 
eling of the viscous terms used in this study and the 
increased streamwise accuracy of the explicit calcu- 
lation because of the large number of streamwise sta- 
tions required. The previous study of Newsome et al. 
used 220 streamwise planes, while the explicit calcu- 
lation made here used 20 220 streamwise planes. The 
same number of points in the vertical direction was 
used in both studies. 
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Figure 4.3. Comparison of computed temperature profiles on a flat plate. = 2.0; Rex, = 1.65 x 10 6 ; x = 0.93. 
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Figure 4.5. Flow field for a hypersonic compression corner. 



Figure 4.6. Comparison of computed pressure coefficients with experimental data. = 14.1; R e L — 1.04 x 10 5 . 
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Figure 4.8. Comparison of computed skin friction coefficients with experimental data. = 14.1; Re/, = 1.04 x 10 5 . 



32 




Normal distance, y 3 Normal distance, 


.4 


.3 




•2 


.1 — 



oL_ 

1.25 


1 .50 

Axial distance, x 


1.75 


2.00 


ire 4.9. Mach number contours on the 15° ramp from explicit upwind solution. M x = 14.1; Re L = 1.04 x 10 4 . 
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Figure 4.15. Hypersonic inlet pressure contour from explicit upwind solution. = 15; Re/, = 8.0 x 10 4 . 
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Figure 4.16. Hypersonic inlet Mach number contour from explicit upwind solution. M = 15; Re/, = 8.0 x 10 4 . 
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5. Three-Dimensional Results 

Two three-dimensional hypersonic viscous flow 
cases were simulated by solving the three-dimensional 
parabolized Navier-Stokes (PNS) equations. The so- 
lution of the three-dimensional PNS equations re- 
quires the approximation of two flux derivatives in 
the crossflow plane. The two-dimensional form of 
the explicit upwind algorithm for the PNS equations 
was extended to solve the three-dimensional PNS 
equations by approximating the convection deriva- 
tives with (1) upwind flux approximations (fully up- 
wind method) and (2) a combination of upwind and 
MacCormack’s (1969) flux approximations (upwind- 
biased method). Both forms of the algorithm are 
used in the first three-dimensional test case for sim- 
ulating hypersonic viscous flow over a cone at high 
angle of attack. A limitation in using the fully up- 
wind method at high angles of attack is discussed. 
The fully upwind form of the algorithm is used in 
the last test case to simulate a Mach 24.5 flow field 
about a generic airplane configuration. A special pro- 
cedure is used for defining the numerical grid in the 
crossflow plane during the development of the sharp 
delta wings. 

The computer code for solving the three- 
dimensional PNS equations used to generate these 
results can be characterized as a MacCormack code 
plus a subroutine to execute the equations in the 
body and appendix of this study for determining 
the flux-difference splitting (FDS). All inner do loops 
were vectorized for use on the NAS Cray 2 computer. 
While the code is highly vectorized, it was developed 
to verify the algorithms rather than to optimize ex- 
ecution speed. The execution times on the Cray 2 
given here can be improved by combining some of 
the calculations and extending the effective lengths 
for multiple do loops. 

5.1. Hypersonic Flow Over a Cone 

The fourth test case simulates laminar, three- 
dimensional hypersonic flow over a 10° half-angle 
cone at an angle of attack a of 24° (fig. 5.1). The 
conditions are 

Moo = 7.95 Re = 4.101 x 10 6 /m L = 0.3048 m 

T tot ai,oo = 755.4 K T w = 309.8 K Pr = 0.72 

a = 0.75 

The above conditions are for the largest angle of at- 
tack considered in Tracy’s (1963) experimental inves- 
tigation. The high angle of attack and free-stream 
Mach number result in a complex flow field because 
of the interaction of the supersonic crossflow with the 
boundary layer. The inviscid flow structure is dom- 
inated by a conical outer shock. The conical outer 


shock’s shape changes on the leeward side of the cone 
because of the growth of the viscous layer. The cross- 
flow is stagnated on the windward side and rapidly 
expands to supersonic speeds as it wraps around the 
circumference of the cone. The boundary layer grad- 
ually thickens as the crossflow moves across the cone 
toward the leeward side. A crossflow separation oc- 
curs as the flow approaches the top of the leeward 
side. The separation region generates an increased 
displacement thickness on the leeward side of the 
cone which expands the position of the outer conical 
shock on the leeward side. A lambda shock forms in- 
side the outer shock wave as the crossflow approaches 
the leeward side to provide the necessary transition 
to subsonic speeds. The complicated flow field is an 
excellent and demanding test case for establishing the 
capabilities of numerical codes. 

The computational grid for the crossflow plane 
consisted of 50 points in the normal direction and 56 
circumferential points. The computational grid for 
the crossflow plane is shown in figure 5.2 using every 
fifth point in the direction normal to the surface. 
The computational grids and solution contours are 
displayed in conical coordinates: 


6y = — arctan 
R 


0z 


Z 

— — arctan 
R 



(5.1) 


where R = yjy 2 + z 2 . The circumferential rays were 
equally spaced around the cone while the grid was 
initially stretched in the normal direction with a 
stretching parameter of 1.12. The outer boundary 
was set outside the expected shock position. The 
grid was made to grow in the marching direction in 
a conical fashion. The code was started from free- 
stream conditions at x = 0.015 with the marching 
step size determined by equation (3.32). The axial 
Mach number varied between 0.2 and 0.5 for the first 
node off the surface at the crossflow stagnation point. 
Typical values of the lawTof-the-wall coordinate (y+) 
for the first point off the wall at solution station 
(x = 0.266) varied from 0.7 to a maximum of 5.0 
at the crossflow stagnation point. 

One reason for using upwind differencing for the 
convection and pressure terms is to eliminate the 
need for additional smoothing or damping parame- 
ters to maintain numerical stability when capturing 
shock waves. For three-dimensional flows solved with 
upwind differences, a problem was encountered hav- 
ing to do with the mapping of a shock wave onto 
a numerical grid. The perfect grid would be ori- 
ented so that the shock wave was contained along 
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an axis in the computational plane. When a shock 
wave does not move uniformly in the computational 
plane, an unrealistic set of initial states occurs when 
the Riemann problem (RP) is solved tangential to 
the shock wave. This results in a nonphysical flux 
error which is not damped out in regions where the 
numerical scheme has low dissipation. For three- 
dimensional flow about a cone, this situation oc- 
curs along the crossflow symmetry line at high an- 
gles of attack. Nonphysical solutions were obtained 
around the crossflow symmetry line for high angles 
of attack when using upwind differencing of the pres- 
sure and convection terms. Low-angle-of-attack cases 
were calculated using upwind differences without the 
above problem. The nonphysical solutions begin to 
develop when the bow shock wave has moved out- 
side the viscous region, usually at 20000 to 30 000 
marching steps from the initial data crossplane. The 
pressure slowly becomes either extremely high or low 
when compared with both the experimental results 
and the next grid point in the tangential direction. 
This problem was solved by two different methods: 
by adding dissipation with Harten’s (1983) device 
when the shocks were not aligned in the computa- 
tional plane, or by differencing the pressure and con- 
vection terms with MacCormack’s (1969) method in 
the direction tangent to a shock wave and using up- 
wind differencing in the normal direction (upwind- 
biased method). The cross section pressure contour 
for the fully upwind method (fig. 5.3(a)) shows the 
difficulty most strongly at the windward crossflow 
symmetry line, while the contour for the upwind- 
biased method (fig. 5.3(b)) shows the correct result 
at the crossflow symmetry. Application of Harten’s 
device was difficult. Harten’s device adds dissipation, 
or smoothing, by artificially preventing the absolute 
value of eigenvalues from decreasing below a certain 
level. To correct the problem at the crossflow sym- 
metry line, the amount of smoothing added resulted 
in a large increase in the boundary layer thickness 
and smearing of the shock wave. The results obtained 
with Harten’s device are not included since they were 
poor compared with the solution obtained with the 
upwind-biased method. The correct solution to the 
three-dimensional flow problem was obtained with- 
out additional damping terms only if a central dif- 
ferencing algorithm (MacCormack’s) in the direction 
tangential to the shock was used. 

The surface pressure distribution at x = 0.325 
(fig. 5.4) agrees well with experimental data on the 
leeward side and is slightly lower than the experimen- 
tal results on the windward side. This is typical of 
previous numerical studies (Gielda and McRae 1986; 
McRae 1976) and has been attributed to the experi- 
mental error associated with the size of the pressure 


tap compared with the boundary layer thickness on 
the windward side. The experimental results include 
cross section surveys of the flow field showing the 
location of the shock wave, viscous region, and min- 
imum pitot tube pressures. These data were taken 
along surface normals at x — 0.2831 (that is, 8.8 cm) 
from the apex of the cone measured along the cone’s 
surface. The numerical results are in a plane nor- 
mal to the centerline of the cone. The numerical re- 
sults were compared with the experimental data by 
projecting the numerical results into a conical coor- 
dinate system. The numerical results for the cross- 
flow plane that bisects the conical experimental data 
are shown in figure 5.5 for the Mach number contour 
at x = 0.266 with the experimental determination 
of shock location, viscous boundary, and minimum 
pitot pressure. Since the flow field is nearly conical 
at this point, the locations of the shock and viscous 
region agree fairly well except on top of the leeward 
side. 

The execution of FDS required 56 percent of the 
CPU time. The new algorithm achieved a compu- 
tational rate of 0.434 x 10 -4 sec per point for one 
complete step of the algorithm. 

5.2. Hypersonic Flow Past a Generic 

Vehicle 

The last test case simulates a laminar hypersonic 
flow field about a generic airplane configuration. The 
purpose of this test case is to demonstrate the capa- 
bility of the upwind algorithm for solving the PNS 
equations to simulate a hypersonic flow field about 
a realistic geometry. The airplane configuration and 
flow conditions are taken from the numerical study 
done by Richardson and Morrison (1987). Note that 
this is a demonstration case since real gas effects are 
not taken into account. The forward part of the 
body is a 4.6° half-angle sharp circular cone, which 
extends 756 in. (19.2 m) from the nose. The cone 
is connected to a cylindrical body which extends to 
1371 in. (34.8 m) from the nose (fig. 5.6). The 12° 
delta wing has a cross section defined by an angle of 
9.327°, set at an angle of attack of 1° relative to the 
fuselage centerline. The delta wing begins 584.6 in. 
(14.8 m) from the nose. The configuration geome- 
try is defined by the following equations where the 
coordinate positions are defined in figure 5.6(b): 

. / 61 in. \ 

fl = mm l I 75T^' 61 m J 

y 0 = —24.5837 in. 

ybw = Vo-(x- 405 in.) tan(l°) - z te 
Vie = Vbw z le = (x- 405 in.) tan(12°) 
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The flow conditions are 

Moo = 24.5 Re = 12 000/in. (4.7 x 10 6 /m) 

L = 1371 in. (34.8 m) T w = 2470°R (1372 K) 

= 490°R (272 K) Pr = 0.72 

a = 1° 

For high Mach number flows at large Reynolds 
numbers, the viscous region becomes extremely thin; 
thus calculation of flow fields around realistic geome- 
tries becomes difficult. The airplane body has a num- 
ber of discontinuous changes in the surface geome- 
try, which have made calculations of flow fields with 
traditional numerical methods difficult. Numerical 
stability problems are often encountered at these lo- 
cations. In the previous study by Richardson and 
Morrison, the thin-layer Navier-Stokes (NS) equa- 
tions were solved with an implicit finite-volume up- 
wind scheme to simulate the flow field. 

The NS solution (Richardson and Morrison 1987) 
was calculated at 52 streamwise stations, with the 
numerical grid at each cross section containing 65 x 65 
points. The calculation of the PNS solution with an 
explicit space marching scheme requires the defini- 
tion of approximately 50000 to 200000 streamwise 
stations. At first, solution to this problem was at- 
tempted with the same numerical grid as the previ- 
ous study, by interpolating (the cross section grids) 
between the 52 streamwise stations to obtain the in- 
termediate stations. This grid development proce- 
dure was satisfactory for space marching the PNS 
equations until the apex of the delta wings was en- 
countered on the fuselage. The point distribution 
and the orthogonality of the numerical grid at the 
body surface required severe changes in the numer- 
ical grid in the streamwise direction at the root of 
the delta wing. The grid points to be used in defin- 
ing the wing surface were collected at the stream- 
wise station immediately before the apex of the delta 
wing. All these points were distributed on the wing 
surface at the next station downstream. Moving all 
the points onto the wing at one time and requiring 
the interior domain to be orthogonal at the surface 
required a large movement of the numerical grid in 
physical space between these two stations. The se- 
vere streamwise change in the numerical grid around 


the wing region was difficult for the space marching 
scheme to handle with realistic space marching step 
sizes. Note that the large streamwise spacing used in 
the NS calculation effectively smooths out any sharp 
changes in the streamwise geometry and the develop- 
ment of the wing. Because of the above problems, a 
different gridding procedure was developed to handle 
the development of the wing by smoothly adjusting 
the physical movement of the numerical grid. 

The new crossflow numerical grids used for this 
case were formed with algebraic stretching functions 
(eq. (4.1)) to cluster the grid at the body surface 
(fig. 5.7). No attempt was made to make the nu- 
merical grid orthogonal at the surface. The point 
distribution on the outer computational boundary 
and the body surface was controlled in an attempt 
to minimize the streamwise changes in the numerical 
grid. The outer computational boundary was a 7° 
cone until the wing had grown large enough to com- 
press the numerical grid to 10 percent of the height 
of the grid along the symmetry line. After this point 
had been reached, the outer boundary was based on 
the linear combination of two cones: a 7° cone at 
the symmetry plane and a cone large enough to in- 
clude the leading edge of the wing plus 10 percent 
of the computational domain at the symmetry line. 
The point distribution on the outer boundary was 
stretched circumferentially so that the points were 
clustered along the plane bisecting the wing’s lead- 
ing edge. The body point distribution before the 
wing is divided into two regions. A constant angular 
point spacing is used above and below the location 
where the apex of the delta wing eventually appears. 
Where the delta wing eventually emerges from the 
fuselage is defined by three points spaced 0.05° apart 
(fig. 5.7). These three points are used to define the 
leading edge of the delta wing. The initial point spac- 
ing on the top and bottom of the wing is equal to the 
point spacing at the leading edge of the wing. As the 
wing grows in size, the points on the body surface are 
rotated one at a time onto the wing in a continuous 
fashion. The amount of rotation is controlled by a 
ratio of the cross section body and wing perimeters. 
When 16 points have been rotated onto both the top 
and the bottom of the wing, the point spacing on 
the wing begins to increase and the rotation stops. 
The point spacing for the three points describing the 
leading edge of the delta wing remains fixed for the 
complete length of the wing. 

The computational grid for half of a crossflow 
plane was defined by 45 points in the 77-direction 
(away from the body surface) and 63 points in £- 
direction (circumferential). Three crossflow grids are 
shown in figure 5.7. The plots of the numerical grid 
show every fifth line in the radial direction. The 
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symmetry boundary condition required 4 of the 63 
points in the (-direction. At the final station, the 
distribution of points in the (-direction was as fol- 
lows: 16 points on the upper fuselage, 12 points 
on the lower fuselage, and 32 points on the wing. 
The outer computational boundary at the symme- 
try plane was equal to the cross section of a 7° cir- 
cular cone. The point distribution along the outer 
computational boundary was algebraically stretched 
(f3 = 1.3) to cluster points around the wing tip. The 
stretching coefficient for the interior domain between 
the body surface and the outer computational bound- 
ary was adjusted to maintain the axial Mach num- 
ber between 0.50 and 0.85 for the first point off the 
surface on the windward symmetry line. The code 
was started using free-stream conditions at 68.55 in. 
(1.74 m) from the apex of the cone. The step size 
was calculated from equation (3.32) using a safety 
factor of 0.95 (inviscid Courant number of approx- 
imately 1.9) up to 585 in. (14.9 m) downstream of 
the cone apex. At this point, the wings began to 
develop and the safety factor had to be lowered to 
0.20 to account for the skewing of the numerical grid 
and the attachment of a shock wave to the leading 
edge of the wing. Once the wings began to develop, 
the pressure boundary condition (eq. (3.31)) had to 
be used to stabilize the surface pressure on the lead- 
ing edge of the wings. The total CPU time on the 
NAS Cray 2 was over 3 hours. It took approximately 
1 hour to advance the solution to the point where 
the wings start, 1 hour to reach the cone-cylinder 
junction, and 1 hour to reach the end of the cylindri- 
cal section. This compares well with the execution 
time for the NS solution of Richardson and Morrison 
(1987), which took approximately 21 hours on the 
Control Data Corp. VPS-32 supercomputer at NASA 
Langley. The code they used executes on the VPS-32 
at roughly the same speed as on the Cray 2. How- 
ever, had the normal grid spacing at the wall been 
as refined as in the Richardson and Morrison com- 
putation, the total execution times would have been 
similar. 

A numerical difficulty was encountered once the 
leading edge of the wing intersected the bow shock 
and moved out into the free stream (at approximately 
700 in. (17.8 m) from the cone apex). A slight pres- 
sure undershoot occurs for the shock captured on the 
ray of points emerging from the leading edge of the 
wing. The pressure undershoot eventually causes a 
numerical instability. The pressure undershoot was 
found to be eliminated by either lowering the value 
of (3 in the flux limiter or by adding explicit smooth- 
ing to the ray of points emerging from the leading 
edge of the wing. The smoothing terms could then 
be removed or (3 increased, but eventually the pres- 


sure undershoot would occur again. The undershoot 
is thought to be caused from the misalignment of the 
grid around the shock wave adjacent to the leading 
edge of the wing. The grid movement causes shock 
wave misalignment with the grid and occurs because 
of expansion of the outer computational boundary 
and adjustment of the stretching coefficient to main- 
tain the streamwise Mach number. A large enough 
grid movement changes the shock position relative 
to the grid points in the computational plane. The 
pressure undershoot is considered to be a result of 
this shock wave position change on the grid. The so- 
lution presented here was calculated with an explicit 
second-order smoothing term added only to the ray 
emerging from the leading edge of the wing. The ex- 
plicit smoothing term added to the solution at these 
points is defined as 

a l(Ei)7/7/ + a 2(Ej)^ (5.2) 

where 

Gl = c lPrm = C2P^ 

The second partial derivatives of (E *) and (p) with 
respect to rj and ( were approximated using second- 
order central differencing. The explicit smoothing 
was applied after the 700-in. (17.8-m) station using 
a coefficient of ci = C 2 = 0.002. 

The PNS results calculated in this study were 
compared with Richardson and Morrison’s (1987) NS 
results. To compare the PNS finite-difference solu- 
tion directly with the NS finite- volume results, the 
cell center locations had to be calculated for the NS 
solutions. Pressure, temperature, and axial velocity 
profiles at the windward symmetry plane are com- 
pared at three stations in figures 5.8, 5.9, and 5.10. 
The profiles agree favorably with respect to values 
before and after the shock, surface values, and pro- 
file shape. Slight differences in the various profiles 
can be partly attributed to the different point spac- 
ing used in the two calculations, the different nu- 
merical integration methods used, and the different 
equation set solved. The first station (256 in. (6.5 m), 
part (a) of each figure) is on the forward part of the 
cone, before the wings appear. The pressure, tem- 
perature, and velocity profiles compare favorably be- 
tween the two methods. The second station (767 in. 
(19.5 m), part (b) of each figure) is 11 in. (0.3 m) 
downstream of the cone-cylinder junction. An ex- 
pansion wave begins to propagate into the flow field 
away from this junction to expand the flow around 
the cone-cylinder corner. The comparison of pressure 
profiles (fig. 5.8(b)) shows that in the NS solution, 
the expansion has propagated farther away from the 
body and is more rounded than the PNS solution. 
The rounding of the expansion wave can be partially 
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attributed to the upstream influence permitted by 
the NS equations. The PNS solution does not round 
the wave since it does not feel the effect of the cor- 
ner until it is reached. Also, the large streamwise 
spacing used in the NS calculation, approximately 
27 in. (0.7 m), effectively smooths out the beginning 
of the expansion wave. The PNS solution predicts 
a slightly higher maximum temperature (fig. 5.9(b)), 
while good agreement is observed between the axial 
velocity distributions (fig. 5.10(b)) at the second sta- 
tion. The third station (1304 in. (33.1 m), part (c) 
of each figure) is near the end of the airplane. The 
strength of the outer bow shock was weakened by 
the expansion of the flow onto the cylindrical body. 
Slightly different shock locations are predicted by the 
two calculations while the surface pressure is in agree- 
ment (fig. 5.8(c)). The PNS solution predicts a small 
drop in pressure outside the edge of the boundary 
layer which is not predicted in the NS solution. This 
drop in pressure could be caused by the turning of the 
flow in the crossflow plane near the saddle point on 
the symmetry plane, as discussed subsequently. The 
numerical grid used in the PNS calculations contains 
almost twice as many points on the lower part of the 
fuselage surface as the numerical grid used in the NS 
calculations. The resolution of the numerical grid 
used in the NS calculation around the lower body 
may not have been adequate to resolve this feature. 
The PNS solution consistently predicts a higher max- 
imum temperature in the boundary layer (fig. 5.9(c)). 

The surface pressures on the wing are compared 
in figure 5.11. A slight difference between solutions 
is observed toward the middle of the wing. A portion 
of this difference may be attributed to the more 
accurate pressure boundary condition used in this 
study or the difference in modeling the leading edge 
of the wing. A solution point is located along the 
leading edge in the PNS calculation while in the NS 
calculation the leading edge is between two cells. The 
pressure contours are compared in figures 5.12, 5.13, 
and 5.14, using 16 identical contour levels for each 
station. The PNS solution more sharply captures 
the bow shock due to the more refined grid at the 
shock location. The PNS prediction of the bow shock 
location agrees with the NS solution. The outer bow 
shock and the decrease in the pressure levels around 
the cone for flow at an angle of attack are clearly 
defined in the pressure contour at station 256 in. 
(6.5 m) for the PNS solution (fig. 5.12). The bow 
shock, expansion wave, delta-wing shock, and change 
in pressure around the cone circumference are shown 
for the PNS solution at station 767 in. (19.5 m) in 
figure 5.13. The pressure contour at station 1304 in. 
(33.1 m) (fig. 5.14) shows the pressure increasing 
from the body to the bow shock and the high pressure 


region created above and below the wing immediately 
inboard the corner shock. 

The details of the flow field in the crossflow plane 
can be investigated using projections of the velocity 
vectors in the appropriate coordinate plane. For lo- 
cations close to the body, the Cartesian plane normal 
to the body can be used to visualize the flow field in 
the crossflow plane. The Cartesian crossflow veloc- 
ity vectors are shown in figure 5.15. The flow field 
next to the lower surface of the fuselage is shown in 
figure 5.16. Note the high crossflow velocity around 
the corner of the fuselage- wing junction, the vortex 
located underneath the fuselage- wing junction, and 
the saddle point located on the symmetry line. The 
pressure drop outside the edge of the boundary layer 
in the PNS solution (fig. 5.8(c)) could be caused by 
the turning of the flow near this saddle point. The 
Cartesian crossflow velocity vectors for the middle 
of the delta wing (fig. 5.15) suggest a strong reverse 
flow. 

The Cartesian plane projection of the velocity 
vectors on the delta wing is misleading since the in- 
viscid flow field on an isolated delta wing is conical. 
To view the flow field on the delta wing, a conical co- 
ordinate system centered at x = 658.9 in. (16.7 m), 
V = —29.02 in. (0.7 m), z = 53.96 in. (1.4 m) on the 
leading edge of the delta wing was used to project 
the velocity vectors onto a conical plane. The loca- 
tion of the conical coordinate system was obtained by 
extending a line along the lower wing- fuselage junc- 
tion until it intersected the leading edge. The use 
of a conical plane to observe the crossflow velocity 
vectors on the wing is an attempt to examine the 
crossflow in a more natural plane. Note that differ- 
ent locations of the conical coordinate system yield 
slightly different results. The conical coordinate sys- 
tem defined above was selected because it represents 
a projection plane that is normal to the leading edge 
and the fuselage-wing junction. The conical veloc- 
ity vectors along the leading edge are shown in fig- 
ure 5.17. The leading edge shock wave and the flow 
separation on the upper and lower surfaces are ap- 
parent from the velocity profiles. The flow separates 
at approximately 9 Z & 7.2° on the top surface and at 
0 Z « 6.2° on the lower surface. The conical velocity 
vectors along the middle of the wing span are shown 
in figure 5.18. Note that the scaling factor for the 
magnitude of the velocity vectors has doubled from 
figure 5.17. The separated flow on the lower surface 
has reattached at approximately 9 Z « 4.6°. The flow 
on the upper surface is more complicated. The flow 
is reattached to the upper surface at 9 Z « 5.8°, and 
the flow separates a second time at approximately 
~ 5.5° with reattachment at 9 Z ^4.5°. 
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The shock waves in the crossflow plane consist 
of a bow shock, leading edge shock, and a corner 
shock. The corner shock is similar to the shock wave 
structure obtained for supersonic flow in the corner 
of intersecting wedges (West and Korkegi 1972). The 
bow shock and the leading edge shock of the wing 
do not intersect, but are joined by a third corner 
shock. The compression corner formed by the cone- 
delta-wing junction is similar to the corner made by 
two wedges. The expansion around the cone-cylinder 
junction weakens the bow shock and complicates the 
flow field around the corner shock. The conical cross- 
flow Mach number contours (fig. 5.19) are parallel 
with the corner shock and indicate that the corner 
shock is essentially conical. The pressure contours 
(fig. 5.20) and the density contours (fig. 5.21) show 
the complicated corner shock wave structure in the 


crossflow plane. A crossflow shock wave is shown 
by the collection of crossflow Mach number contours, 
pressure contours, and density contours on the upper 
wing surface centered at approximately 9 y ~ —0.5°, 
6 Z « 5.5°. The interaction of the expansion wave 
has diffused the effect of the refracted leading edge 
and bow shock wave shown in the pressure contour 
at station 767 in. (19.5 m) (fig. 5.13). 

To demonstrate the complete flow field more 
clearly, the pressure and Mach number contours are 
shown in color in figures 5.22 and 5.23. The sepa- 
rated flow on the wing and the vortex underneath the 
wing-fuselage junction can easily be seen in the color 
Mach number contour. The high pressure region in- 
board of the corner shock and the gradual pressure 
decay to the body surface from the outer shock wave 
are shown in the pressure contour. 



Shock wave 



Figure 5.1. Hypersonic flow over a cone at an angle of attack. 
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Figure 5.4. Comparison of computed surface pressures with experimental data. = 7.95; a = 24°; Re x = 4.2 x 10 5 
7V T total,oo = °- 4 0; x = °- 325 - 



Figure 5.5. Comparison of computed Mach number contours with Tracy’s flow field survey. M 0 0 ~ 7.95; a = 24°; 
Re x = 3.6 x 10 5 ; x = 0.266. 
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Figure 5.6. Generic airplane configuration. 
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(b) x = 767 in. (19.5 m). 
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(b) x = 1304 in. (33.1 m). 


Figure 5.7. Computational grids for three crossflow planes of the airplane. 
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(c) x = 1304 in. (33.1 m). 

Figure 5.8. Comparison of computed pressure profiles on the windward symmetry plane. = 24.5; Re = 12 000/in. 
(4.7 x 10 6 /m); a = 1°. (1 in. = 0.0254 m.) 
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Figure 5.9. Comparison of computed temperature profiles on the windward symmetry plane. Moo 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. (1 in. = 0.0254 m.) 
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Figure 5.10. Comparison of computed axial velocity profiles on the windward symmetry plane. = 24.5; 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. (1 in. = 0.0254 m.) 
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Figure 5.11. Comparison of computed surface pressures on the wing at station 1304 in. (33.1 m). M x — 24. 
Re = 12 000/in. (4.7 x 10 6 /m); a = 1°. (1 in. = 0.0254 m.) 
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Figure 5.12. Comparison of computed pressure contours at station 256 in. (6.5 m). M (*> = 24.5; 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. 
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(b) PNS explicit upwind solution. 

Figure 5.14. Comparison of computed pressure contours at station 1304 in. (33.1 m). M» = 24.5; 
Re = 12 000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.16. Cartesian crossflow velocity vectors on the windward side of the fuselage at station 1304 in. (33.1 m). 
A/oo = 24.5; Re = 12000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.17. Conical crossflow velocity vectors near the leading edge of the delta wing at station 1304 in. (33.1 m). 
Moo = 24.5; Re = 12 000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.18. Conical crossflow velocity vectors on the delta wing at station 1304 in. (33.1 m). M 0 0 = 24.5; 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.19. Conical crossflow Mach number contours at station 1304 in. (33.1 m). M x — 24.5; 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.20. Computed wing pressure contours at station 1304 in. (33.1 m). M m = 24.5; 
Re = 12 000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.21. Computed wing density contours at station 1304 in. (33.1 m). M x = 24.5; Re = 12000/in. (4.7 x 10 6 /m); 
q = 1°. 


59 



Nondimens ional 
pressure 


i 

l 




0.00104 

.00146 

.00189 

.00231 

00274 

.00316 

.00358 

.00401 

.00443 

00486 

.00528 

.00571 

00613 

00656 

00698 

.00740 

00783 

00825 



Figure 5.22. Computed pressure contours at station 1304 in. (33.1 m). M c c = 24.5; 
Re = 12 000/in. (4.7 x 10 6 /m); a = 1°. 
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Figure 5.23. Computed Mach number contours at station 1304 in. (33.1 m). M 0 c = 24.5 
Re = 12000/in. (4.7 x 10 6 /m); a = 1°. 
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6. Concluding Remarks 

A new algorithm for solving the three-dimensional 
parabolized Navier-Stokes (PNS) equations has been 
developed. The new algorithm is an explicit finite- 
difference scheme which uses upwind flux approxi- 
mations for the pressure and convection terms and 
central differencing for the viscous and heat flux 
derivatives. The upwind flux approximations for 
the pressure and convection terms are based on 
the solution of an approximate Riemann problem 
(RP) for the PNS equations using a modification of 
the method proposed by Roe for steady supersonic 
flow of an ideal gas. Roe’s method is extended to 
solve an approximate RP in E* space for the three- 
dimensional PNS equations transformed into gener- 
alized coordinates and to include the subsonic pres- 
sure splitting technique of Vigneron. The algorithm 
is shown to capture strong shock waves without addi- 
tional damping terms that depend on adjustment of 
solution-dependent coefficients. The execution time 
for the new algorithm is approximately the same as 
a central difference code, since the upwind differenc- 
ing of the pressure and convection terms uses ap- 
proximately 50 percent of the central processing unit 
time and doubles the Courant-Friedrichs-Lewy sta- 
bility limit. The algorithm has proven to be efficient 
for use on vectorized computing machines since all 
inner and some outer do loops are vectorized. 

The new algorithm is demonstrated for two- and 
three-dimensional supersonic and hypersonic lami- 
nar flow test cases. The test cases agree favorably 
with both experimental data and numerical results 
obtained using other numerical methods. Accurate 
flat plate boundary layer profiles are calculated us- 
ing the new algorithm started from free-stream initial 
conditions. Previous numerical calculations by oth- 
ers of hypersonic flow over a ramp have demonstrated 
nonphysical pressure oscillations in the solution. The 
new algorithm clearly resolves the pressure field for 
the ramp without this difficulty. A complicated in- 
let flow field containing intersecting and reflecting 
shock waves is computed and demonstrates the ro- 
bust shock capturing obtained with upwind flux ap- 
proximations of the pressure and convection terms. 


Multidimensional numerical algorithms using up- 
wind differencing of the pressure and convection 
terms have been applied with additional dissipation 
terms to eliminate a loss in accuracy at certain loca- 
tions in the flow field. The flow field over a cone at 
an angle of attack of 24° is computed using a combi- 
nation of upwind and MacCormack differencing for 
the pressure and convection terms. This combination 
scheme eliminated a loss of accuracy at the symme- 
try plane without additional numerical dissipation. 
The flow field about a generic hypersonic airplane at 
Mach 24.5 and an angle of attack of 1° is calculated 
using the new algorithm. This flow field has not been 
previously solved using a noniterative space marching 
method. A special algebraic grid generation routine 
is used which eliminated difficulties associated with 
the numerical grid at the apex of the delta wings. 
In summary, the numerical results obtained with the 
new algorithm more clearly and accurately resolve 
the flow field features than previous results obtained 
with other methods for solving the PNS equations. 

The research performed in the course of this study 
produced the following additional significant results: 

1. The eigenvectors were determined for use in 
solving the approximate RP in E*-space for 
the three-dimensional PNS equations trans- 
formed into a generalized coordinate system 
and including Vigneron’s splitting of the.sub- 
sonic pressure gradient. 

2. The difficulty associated with applying Roe’s 
method in the subsonic region with a non- 
iterative space marching scheme for solving 
the PNS equations including Vigneron’s pres- 
sure splitting procedure was identified and 
overcome. 

3. A simple method was developed for modifying 
the one-sided differencing in MacCormack’s 
method into an upwind differencing scheme. 

4. An increase in the stability of the scheme 
was obtained by solving for the value of 
Vigneron’s pressure splitting coefficient using 
a cubic equation in terms of the dependent 
flux vector. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
November 15, 1990 
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Appendix A 


the above transformation are 


Generalized Transformation 


A general transformation was used in this study 
to transform the governing equations from the phys- 
ical domain (:r,y,z) to the computational domain 
(£, V) 0* The transformation is made so that the gov- 
erning equations can be solved on a uniformly spaced 
computational grid. One of the advantages of gener- 
alized transformation is that it eliminates the need 
to interpolate the body surface onto the numerical 
grid. The transformation is of the following form: 


Zx = 

£y = £z = 0 

Vx = -J{vzz c _ vc z d 

T)y = Jx £Z C 
Vz - 

Cx = J{y^Zr, - yr,Z c ) 

C y ~ Zjj 

Cz = Jx^Vrj J 


(A.2) 


where the subscripts indicate differentiation and J is 
the Jacobian of the transformation 

J = y (A.3) 

H \Vri z Q - y<: z v) 


* = m ' 

T] = T]{x,y,z) > 
C = C (x,y,z) . 


(A.l) 


The derivations of the formulas for the metrics 
(£x, r)x ) Vyi Vz , Cx, C y, Cz) of a generalized transfor- 
mation are given in the text by Anderson, Tannehill, 
and Pletcher (1984). The formulas for the metrics of 


Gielda and McRae (1986) have shown that for 
the PNS equations solved by MacCormack’s (1969) 
method, the geometric conservation law (GCL) terms 
are not zero , for any combination of possible differ- 
encing of the metrics. Therefore, the metrics are cal- 
culated once for each new space marching step us- 
ing a single differencing approximation. The partial 
derivatives of x, t/, and 2 with respect to 7] and C are 
numerically formed using second-order central differ- 
ences. The partial derivative x ^ is approximated with 
a first-order backward difference. 
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Appendix B 

Eigenvalues, Eigenvectors, and Wave Strengths for the Three-Dimensional PNS 
Equations 

The formulas used in constructing the approximate Riemann solution are given in this appendix. The 
formulas are valid for a Riemann problem (RP) constructed in either the £-77 or the plane when the 
appropriate metrics are used. The index j is used to represent the points in the plane in which the RP is 
being solved. The determination of the eigenvalues and the eigenvectors shown here for the three-dimensional 
inviscid PNS equations in generalized coordinates was accomplished in part using the symbolic manipulation 
language MACSYMA. 


B.l. Square- Root Averaging 

The matrix A defined by equation (3.7) is formed from square- root- averaged variables so that conservative 
properties are maintained. The square- root-averaged variables result from averaging a special parameter vector. 
The parameter vector p has properties which not only maintain conservation but also simplify the matrix 
algebra used in determining the eigenvalues and eigenvectors (see Roe 1981). The averaged components of p 
are defined as the square-root-averaged Cartesian velocities and enthalpy. The definition of p is 


Pj = ypj [*> v 'r w j> h J 


iT 


(B.l) 


where p', u', v\ w\ h! are the original unaveraged density, Cartesian velocities, and enthalpy. The square- 
root-averaged velocities and enthalpies are 



u . . 1 


v ■ , 1 




0-5(P2j+1 + P2j) 
0-5(pij+i + Pij) 


•^+1+^ 


L- 1 1+1 
J+5 


iw 


’j + 1 


+ Wj 


+ 1 


i h[ 


l j + 1 + h 'j 


4- 1 




*j+^+i + u i 

R i+h + 1 




(B.2) 


where the first subscript of p indicates the vector component. Unless otherwise stated, all the variables that 
are defined in the following are understood to be at j 4 - \ and are formed using the above defined square-root- 
averaged variables. 


B.2. Metrics 

The metrics are defined so that if the nonmodified flux is forward (so the modified flux term contains 
elements based on the positive eigenvalues) and all the eigenvalues are positive, the resulting difference is 
backward. This happens if we use the metrics at j + 1 to form df + and at j to form df”. For example, if the 
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RP is being solved in £-77 plane, the metrics would be defined as 


Vx 


W+i “ J 

m 

. 1 

J+2 J 

m 

= T ^- 


j + h J 

m j 


(B.3) 


where m = j + 1 for df + and m — j for df . For the plane, 77 would be replaced by £ and the differencing 
would be done in the other plane. 


B.3. Eigenvalues 

The eigenvalues for the three-dimensional, inviscid PNS equation set in generalized coordinates are 


Ai 


~ a 2 ~ 


yja\ -4aia 3 
2a\ 


A2 = A3 = A4 


v 

u 


(B.4) 


A5 = 


—CL2 + a 2 — 4aia3 


2a\ 


where 


v = n x u + riyV + n z w 

ai = u>(u 2 — (?) + 7(1 — lj)v? 

a >2 — (a; + 1)(— vu + n x (?) - 7(1 — u>)vu 

03 = V 2 - c 2 (nl + n.y + n 2 z ) 


c 2 = (7 - 1) 


h—\{u 2 + v 2 + w 2 ) 


and u) is Vigneron’s (1978) coefficient and 7 is the ratio of specific heats. 


B.4. Eigenvectors 

The eigenvectors for the three-dimensional, inviscid PNS equations set in generalized coordinates are 
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.-h. 
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h 
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where 


w = n y w — n z v n't = riy + v\ 
q 2 = u 2 + v 2 4 - w 2 vt = v — n x u 


v g = uv — n x u 


C 2 — a\ — 4 aia 3 


Af = (1 - uj) [(7 - 1 )uv + n x c 2 j 

V = 2 jo; [(u 2 - c 2 )n 2 + v 2 j + 7 u(l - lj)(uti 2 - n x vt ) j 


Rn — 


C+Af 

V 


Ftp — 


C-AT 

V 


B.5. Wave Strengths 

The wave strengths can be determined by solving equation (3.11). This results in the following: 

hAEt-AEZ 
a2 = 2 h~^ 

s = n$u AE% + Vg (% AE% + n z A £|) 
rfu 2 + VgVt 


013 = 


Tly A_£/^ 71 2 A Eg — wS 

2 rtfh 


a 4 = 


as w — A El -f a 2 + S 



(B.6) 


a l 


(C-Af)(S- 2q 4 ) V [v t A El - u(riy A + n z AEj)] 
2£ 2£(n^u 2 + VgVt) 


as = 5 — c*i — 2a4 
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